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Abstract. Numerical simulations of two-dimensional capillary filling using the pseudo-potential lattice 
Boltzmann model for multiphase fiuids are presented, with special emphasis on the role of finite-vapour 
density effects. It is shown that whenever the density of the light-phase exceeds about ten percent of the 
dense phase, the front motion proceeds through a combined effect of capillary advection and condensation. 
As a result, under these conditions, the front proceeds at a higher speed as compared to the Washburn 
prediction. It is suggested that such an acceleration effect might be observed in experiments performed 
sufficiently close to critical conditions. 
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1 Introduction 



In the recent years, increasing attention has been paid 
to the use of lattice Boltzmann (LB) techniques beyond 
the strictly hydrodynamic regime, i.e. for micro and nano- 
fiuidic applications. LB holds great potential for these 
applications because it can offer an optimal compromise 
between the physical realism of atomistic methods and 
the computational efficiency of continuum fluid mechan- 
ics. However, turning this potential into a robust and re- 
liable prediction tool for complex micro/nanofluidic prob- 
lems, requires a detailed and quantitative validation pro- 
gram on a number of test and benchmark problems. In 
this work, we present one such validation study, as ap- 
plied to the important problem of capillary filling. Cap- 
illary filling is an old problem, originating with the pi- 
oneering works of Washburn yy and Lucas [2j. Recently, 
with the explosion of theoretical, experimental and numer- 
ical works on microphysics and nanophysics, the problem 
attracted considerable renewed interest [3ji4ji5ii6j . Capil- 
lary filling is a typical "moving contact line" problem, 
in which interface motion can only take place thanks to 
subtle non-hydrodynamic effects occurring at the contact 
point between liquid-gas and solid phase. As a result, it 
offers an ideal testground for assessing the capabilities of 
mesoscopic methods, such as LB, to describe phenomena 
beyond the hydrodynamic regime. 



2 The Washburn description 

A simple and yet powerful picture of capillary motion was 
developed back in 20 's with the pioneering work of Wash- 
burn and Lucas. The Washburn picture is easily derived 
by writing the momentum equation for the mass of moving 
liquid in the capillary. This leads to the following equation: 



2"fcos{9) 12fiizz 



H 



(1) 



where z{t) is centerline position of the advancing front 
at time t, H the transverse size of the channel, 7 is the 
surface tension, 9 is the static contact angle, pi the density 
and ^ii the viscosity of the liquid. The factor 12 stems 
from the specific geometry considered here, i.e. two infinite 
parallel plates separated by a distance H - see fig. [1]). 
More precisely, in the overdamped regime, the Washburn 
solution reads as follows: 



jHms(e)_~ 



(2) 



To the purpose of comparing with experimental results, it 
is useful to recast this relation in dimensionless form (re- 
duced units), t = t/tcap and z = z/ H , being the capillary 
time tcap ~ Hfj,/"f. This leads to the universal law: 
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The Washburn solution ([3]) above spells troubles for mi- 
crofluidic applications where fast-fill is a goal. Indeed, ([3|) 
clearly shows that the fi'ont moves at increasingly lower 
speed as it propagates down the capillary. The reason is 
easily understood by noting that, according to the eq. 
[U the front position z{t) evolves like the trajectory of a 
point-like particle, whose mass grows linearly with the po- 
sition z, driven by a constant capillary force and damped 
by a dissipative force also growing linearly with z. By ne- 
glecting inertia (overdamped regime), the capillary force 
competes with a linearly increasing dissipation, which ul- 
timately leads to a front speed z decaying like 1/z. Strate- 
gies to do away with this basic limitation make the object 
of intense scientific and technological research, such as in- 
terface functionalization, optimal coating, and related is- 
sues. As already remarked in the literature [5j, the asymp- 
totic behaviour ([3]) hinges on a number of simplifying as- 
sumptions, namely: (i) the inertial terms in the Navier- 
Stokes equation are negligible, (ii) the instantaneous bulk 
profile is given by the Poiseuille flow, (iii) the microscopic 
slip mechanism which allows for the interface motion is not 
relevant for bulk quantities (such as the overall position 
of the interface inside the channel), (iv) inlet and outlet 
phenomena can be neglected (limit of infinitely long chan- 
nels); (v) the liquid is filling in a capillary, either empty or 
filled with gas whose total mass is negligible with respect 
to the liquid one. None of these assumptions needs to be 
true in actual experiments, and to the purpose of com- 
parison with experimental data, it is therefore important 
to address these limitations within a suitably generaliza- 
tion of the Washburn equation, to which we shall return 
shortly. In this paper, we shall be interested specifically in 
the finite- vapour density issue (v). 



3 LBE for capillary filling 

The model used in this work is a suitable adaptation of 
the Shan-Chen pseudo-potential LBE [7] with hydropho- 
bic/hydrophilic boundaries conditions, as developed in [8l 
[9j. Other models with different boundary conditions and/or 
non- ideal interactions have been also used in [lOj . 

The geometry is depicted in fig. ([1]). The bottom and 
top surface is coated only in the right half of the channel, 
with a boundary condition imposing a given static contact 
angle [8]; in the left half, we impose periodic boundary 
conditions at top and bottom surfaces, so as to realize a 
flat liquid-gas interface mimicking an "inflnite reservoir" . 
Periodic boundary conditions are also imposed at the two 
lateral sides, so as to ensure total conservation of mass 
inside the system. 



3.1 LBE algorithm for multi-phase flows 

We start from the usual lattice Boltzmann equation with 
a single-time relaxation jllill2| : 

At 

TB 

(4) 
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Fig. 1. Geometrical set-up of the numerical LBE. The 2 di- 
mensional geometry, with length 2L and width is divided in 
two parts. The left part has top and bottom periodic bound- 
ary conditions such as to support a perfectly fiat gas-liquid 
interface, mimicking a "infinite reservoir". In the right half, of 
length L, there is the true capillary: the top and bottom bound- 
ary conditions are those of a solid wall, with a given contact 
angle Q jSj. Periodic boundary conditions are also imposed at 
the west and east sides. 



where fi{x,t) is the kinetic probability density function 
associated with a mesoscopic velocity c;, tb is a mean 
collision time (with At a time lapse), //'^^''(p, pu) the equi- 
librium distribution, corresponding to the Maxwellian dis- 
tribution in the continuum limit. From the kinetic distri- 
butions we can deflne macroscopic density and momentum 
fields as pT|[T2]: 



Pix)^^fi{x); pu{x) =^cifi{x) 



(5) 



For technical details and numerical simulations we shall 
refer to the nine-speed, two-dimensional 2DQ9 model 
The equilibrium distribution in the lattice Boltzmann equa- 
tions is obtained via a low Mach number expansion of the 
equilibrium Maxwellian [lllll2j . In order to study non- 
ideal effects we need to supplement the previous descrip- 
tion with an interparticle forcing. This is done by adding 
a suitable in ([1]) . In the original model [7] , the bulk in- 
terparticle interaction is proportional to a free parameter 
(the ratio of potential to thermal energy) , G, entering the 
equation for the momentum balance: 
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w{\ci\^)^{x,t)i^{x + ciAt,tyi (6) 



being w{\ci\^) the static weights for the standard case of 
2DQ9 TI] and %lj{x,t) = ilj{p{x,t)) the pseudo-potential 
function which describes the fluid-fluid interactions trig- 
gered by inhomogeneities of the density proflle (see [IllHl 
E] for details). 

One may show [8,9J that the above pseudo-potential, 
leads to a non-ideal pressure tensor given by (upon Taylor 
expanding the forcing term): 



P — 



(7) 
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where Cg is the sound speed. This approach allows the def- 
inition of a static contact angle 9, by means of a suitable 
value for the pseudo-potential = i^iPw) [S]j which can 
span the range G [0° : 180°]. Moreover, it also defines 
a specific value for the surface tension, 7/g, via the usual 
integration of the offset between normal and transverse 
components of the pressure tensor along the liquid-gas in- 
terface [7„8^. 

As to the boundary conditions on the Boltzmann popu- 
lations, the standard bounce-back rule is imposed. One 
can show that the bounce-back rule gives no-slip bound- 
ary conditions up to second order in the Knudsen num- 
ber in the hydrodynamical limit of single phase flows |14j . 
In the presence of strong density variations, close to the 
walls across the interface, the velocity parallel to the wall 
may develop a small slip length (of the order of the in- 
terface thickness, oc ^) which, in turn, allows the inter- 
face to move. At variance with continuum formulations, in 
which the slip length is prescribed a-priori, our mesoscopic 
model allows slip flow to develop on its own as a result of 
the fluid-wall interactions coded in the mesoscopic wall 
potential i/'u'- However, it is difficult to control exactly 
this phenomenon, because even imposing an exact no-slip 
boundary condition at the wall [15 , the model does de- 
velop its own dynamics at the flrst node away from the 
wall, thus leading to an overall non-zero slip velocity. To 
the purpose of controlling the capillary flUing, one may 
reabsorb all these effects within the usual Maxwell slip 
boundary conditions: Us = XgdnU. It is easy to show that 
in presence of a slip velocity, the Poiseuille profile becomes: 



u{y) = 6 



u y{H-y) + XsH 
i?2 1 + 6X,/H 



(8) 



where the velocity of the front must be identified with the 
mean velocity, u = 1/H u{y)dy — z. Therefore, the 
Washburn law ^ becomes: 



where 



A = l + 6 



H 



(9) 



(10) 



accounts for finite-slip slow effects. Under ordinary condi- 
tions, the slip length is of the same order of the atomistic 
interaction scale, hence £,/H < 10~^, so that slip-effects 
can be neglected to all practical purposes. This situation 
can drastically change in the presence of superhydrophic 
effects, although we shall not be concerned with these 
problems in the present work. 



3.2 Generalized Washburn equation 

As already remarked many years ago [TB], the Washburn 
law ([3]) holds in the limit where inertial forces can be ne- 
glected with respect to the viscous and capillary ones. This 
cannot be true in the early stage of the filling process, 
where strong acceleration drives the interface inside the 



capillary. However, putting typical numbers for microflu- 
idic devices {H ~ l/im, 7 ~ O.Q72N/m, pi ~ 10~^kg/m^, 
fii ~ 10^^ Ns/m'^), it is readily checked that the transient 
time, Tdiff = H^pi/pf, is of the order of a few nanosec- 
onds, hence completely negligible for most practical pur- 
poses. 

Another important effect which must be kept in mind 
when simulating capillary flUing, is the unavoidable "re- 
sistance" of the gas occupying the capillary during the 
liquid invasion. This is a particular "sensitive" issue for 
LB techniques, because reaching the typical 1 : 1000 den- 
sity ratio between liquid and gas of experimental set up, 
represents a challenge for most numerical methods, par- 
ticularly for multiphase Lattice Boltzmann, typically op- 
erating in the regime 1 : 10 to 1 : 100. In order to take 
in to account both effects, inertia and gas dynamics, one 
may write down the balance between the total momentum 
change inside the capillary and the force (per unit width) 
acting on the liquid+gas system (here untilded symbols 
denote physical units): 



d{iM{t)) 
dt 



cap 



(11) 



where M (t) = Mg + Mi is the total mass of liquid and gas 
inside the capillary at any given time. The two forces in the 
right hand side correspond to the capillary force, Fcap = 
2'^cos{9), and to the viscous drag Fyis = —2{ng{L — z) + 
liiz)dnu{Q). Following the notation of fig. ([I]) and the ex- 
pression for the velocity profile ([S]) one obtains the final 
expression (see also |17) for a similar derivation, without 
considering the slip velocity): 

{pg{L - z)^ piz)z + {pi - pg){zf = 

jcom Uz (12) 

In the above equation for the front dynamics, the terms 
in the LHS take into account the fluid inertia. Being pro- 
portional either to the acceleration or to the squared ve- 
locity, they become negligible for long times. Washburn 
law plus the slip correction ^ is therefore correctly re- 
covered asymptotically, for t — > 00, and in the limit when 
Pgl pi — >■ 0. The above equation is exact, in the case where 
evaporation-condensation effects are negligible, i.e. when 
the gas is pushed out of the capillary without offering 
any mechanical resistance to the advancing liquid. This 
is not the case for most mesoscopic models available in 
the literature [71IT5] , based on a diffusive interface dynam- 
ics [T^. As we shall see, only when either the limit of 
thin interface ^/-ff is reached or when the gas phase 
density is negligible, Pgj pi 0, the dynamics given by 
is correctly recovered. Otherwise, deviations induced 
by condensation/evaporation effects, are observed, which 
may result in significant departure from the Poiseuille pro- 
file inside the gas phase. 



4 Numerical results 

In a previous paper [3T], we pointed out the effects of 
finite-density jumps, interface width and inertia in LB 
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Fig. 2. Evolution of the front coordinate, z{t), for two values 
of the density ratio, pi/pg — 34 and pi/pg — 11. For the sake 
of comparison, the numerical solution of the generalized Wash- 
burn equation with the same value of pi/pg is also reported. 
The figure clearly shows that only the case pi/pg — 34 provides 
good agreement between the two. 

simulations of capillary filling. In the present work, we 
keep focus on the details of the former alone, leaving a 
more complete discussion to a future and lengthier publi- 
cation. 



4.1 Finite-density effects 

In Figure [2] we show the front trajectory z{t) for a given 
resolution H — 121 and two different density ratios, pi/pg = 
11, (G = 5.0), pi/pg = 34, (G = 6.0). The LB solution is 
compared against numerical solution of the generalized 
Washburn equation (fT2|) . From this figure, it is clearly 
appreciated that the Washburn solution is quantitatively 
reproduced only for density ratios above 30. 

It is instructive to inspect the internal structure of the 
flow (Fig. (3]). In the left panel, the expected Poiseuille 
flow is clearly visible for both liquid and gas phases. Of 
course, distortions from such asymptotic behaviour take 
place in the vicinity of the moving interface, but since the 
vapour phase is light enough, such distortions are rapidly 
reabsorbed away from the interface. As the vapour phase 
is made heavier, however, such reabsorption no longer 
takes place, and a qualitative change is observed in the 
flow pattern, with a clear emergence of a 'anti-Poiseuille' 
flow in the gas phase. Such anti-Poiseuille flow appears 
to be a direct consequence of the fact that, within this 
parameter regime, the front evolves under the combined 
effect of two distinct mechanism: capillary drive and front- 
condensation. The extent to which such condensation holds 
in real experiments, remains an open issue at this point, 
and one to which LB can hopefully contribute useful in- 
sights for the future. Here, we only wish to offer a few 
remarks. First, we observe that since front condensation 
proceeds diffusively, it contributes the same t^^'^ scaling 
exponent as capillary advection. This explains why one 
can recover the correct 1/2 Washburn exponent, and yet 
miss the right prefactor. This is relevant to LB, since it is 
known that LB simulations are affected by spurious cur- 




Fig. 3. Streamwise velocity field for the case pi/pg = 11 (left) 
and pt/pg = 34 (right). The dark region marks the dense phase. 

rents due to lack of isotropy of the higher-order kinetic 
moments. 
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Fig. 4. Two-dimensional velocity pattern for the case pi/pg — 
11 (upper fig.)and pi/pg — 34 (lower fig.). The presence of 
two rolls (vortices) near the solid wall, ahead and past the 
front is clearly visible. Equally visible, is in the upper figure, 
the inverted flow in the light phase, far ahead of the front. The 
far-field velocity ahead of the front of the lower figure, however, 
recovers the Poiseuille profile. 

One might argue that the rolls shown in figures [H 
are just a discreteness artifact due to spurious currents. 
We are reasonably confident that this is not the case, be- 
cause the same phenomenon has been observed in inde- 
pendent LB simulations using the momentum-conserving 
free-energy formulation, which is known to suffer less se- 
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vere spurious-current effects |20| . Moreover, by running 
the same case at doubled resolution, the same picture has 
been obtained, namely the vortex intensity was found to 
stay nearly unchanged. Since it is know that spurious cur- 
rents decay with grid resolution this provides further 
confidence that front condensation is a real physical ef- 
fect which may occur in experiments near the liquid-gas 
phase-transition. However, since in LB simulations the in- 
terface width is far thicker than in actual experiments, 
it is possible that the condensation effects predicted by 
the simulation might significantly overestimate the real 
ones. It would therefore be desirable to develop a theoret- 
ical model to predict the share of condensation/capillarity 
as a function of the liquid-gas density ratio, possibly by 
extending the generalized Washburn equation along the 
ideas outlined in |22j. Work along these lines, as well as to 
understand the role of dynamic versus static contact an- 
gle and additional dissipative effects induced by the front 
deformation during propagation, is currently underway. 



5 Conclusions 

The present study shows that Lattice Boltzmann mod- 
els with pseudo-potential energy interactions are capa- 
ble of reproducing the basic features of capillary filling, 
as described within a (generalized) Washburn approxima- 
tion. Two conditions for quantitative agreement have been 
identified in the recent past: i) a sufficiently high density 
contrast between the dense/light phase, pi/pg > 10 and 
a sufficiently thin interface, £^/H < 0.1. In this paper we 
have focused on the former effect, and shown that at suffi- 
ciently high vapour density (above 1/10 of the liquid den- 
sity) front motion proceeds through a concurrent combi- 
nation of capillary advection and front condensation. The 
latter is sustained by an inverted Poiseulle motion in the 
light phase, providing a steady supply of vapour to the 
front location, where condensation takes place. Indepen- 
dent LB simulations, based on a different LB formulation, 
do confirm the same picture, thereby lending weight to 
the hypothesis that at finite vapour density, capillary fill- 
ing may proceed through a combination of surface ten- 
sion drive and condensation. It is also observed that since 
both mechanisms feature the same scaling exponent, such 
an hypothesis can only be tested through a detailed and 
quantitative analysis of the prefactors and not just the ex- 
ponents. Our data suggest that motion by condensation 
proceeds faster than the capillary one. Should this indica- 
tion receive experimental confirmation, one might devise 
ways of accelerating capillary filling by performing exper- 
iments near the critical point. 
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